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\ SUMMARY 

\ introductory  account  is  given  of  the  use  of  isoparametric  elements  infinite  element 
structural  analysis.  The  basic  theory  is  described  and  two  simple  examples  are  worked  out 
in  detail;  these  both  relate  to  the  torsion  of  bars.  Brief  mention  is  also  made  of  other 
applications  of  isoparametric  elements. 
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NOTATION 


an  Term  of  element  matrix 

bi  Constant  term  associated  with  element  matrix 

tn,  V(,  wt  See  Equations  (3.4.5);  also  Tables  2 and  3 
x,  y Co-ordinates  in  physical  plane 

xu  yt  Co-ordinates  of  nodal  points  in  physical  plane 

A,  B,  C See  Equations  (3.4.3) 

G Shear  modulus 

/ Minimisation  integral  for  cross-section 

It  Minimisation  integral  for  element 

J Jacobian  of  mapping  function 

K Torsion  constant  for  cross-section 

Ki  Torsion  constant  for  element 

P(  Interpolation  function  for  four-node  element 

Qt  Interpolation  function  for  eight-node  element 

-q  Co-ordinates  in  transformed  plane 

tx,  Ty  Shear  stresses  in  plane  of  cross-section 

0(f,  rj)  Torsion  function  in  transformed  plane 

’F(x,  y)  Torsion  function  in  physical  plane 

Nodal  values  of  torsion  function 


1.  INTRODUCTION 


Almost  all  aircraft  structural  analysis  these  days  is  carried  out  using  computer  programs 
based  on  the  finite  element  method  of  analysis.  However,  at  least  for  the  novice  user  of  such 
programs,  there  is  sometimes  an  uncertainty  about  the  properties  of  the  various  types  of  elements 
incorporated  in  the  program  libraries.  This  seems  to  apply  particularly  to  the  so-called 
“isoparametric  elements”  even  though  these  were  developed  around  a decade  ago1,2  and  have 
received  some  attention  in  the  text  book  literature.1 2 3 

Perhaps  one  of  the  reasons  why  uncertainties  exist  in  the  case  of  isoparametric  elements  is 
that  even  simple  problems  involving  them  can  lead  to  formidable  calculations  and  in  few  places 
can  one  find  illustrative  examples  worked  in  detail.  It  is  with  these  thoughts  in  mind  that  the 
following  expository  account  of  isoparametric  elements  has  been  prepared.  Here,  after  a review 
of  the  basic  principles,  two  illustrative  examples  are  presented  at  length.  In  order  to  keep  the 
calculations  of  reasonable  size  both  examples  relate  to  the  torsion  of  bars,  this  being  one  of  the 
simplest,  but  still  non-trivial,  problem  types  for  which  isoparametric  elements  can  be  usefully 
employed.  (The  torsion  problem  is,  for  example,  simpler  than  the  plane  stress  problem 
associated  with  sheet  structures.)  Again,  to  keep  the  calculations  sufficiently  concise  so  that  they 
can  be  displayed  at  length,  only  crude  finite  element  subdivisions  are  employed  and  no  pretence 
is  made  about  the  precision  of  the  answers  for  the  illustrative  examples.  Finally,  an  outline  is 
given  of  other  applications  of  isoparametric  elements,  including  their  use  for  determining  the 
stress  intensity  factor  in  fracture  mechanics  problems. 


2.  PRINCIPLES  OF  ISOPARAMETRIC  ELEMENTS 


2.1  General 


The  broad  approach  in  any  finite  element  analysis  is  to  consider  the  particular  region  of 
interest  (which,  in  the  illustrative  examples  to  be  discussed  later,  will  be  the  cross-section  of  a 
bar  in  torsion)  as  subdivided  into  a finite  number  of  smaller  regions,  or  “elements”,  these 
elements  being  deemed  connected  at  certain  points  or  “nodes”.  In  order  to  determine  how  some 
function  of  interest  (which  in  the  illustrative  examples  will  be  the  Prandtl  torsion  function) 
varies  throughout  the  complete  region,  a relatively  simple  form  of  variation  is  assumed  over 
each  of  the  elements.  The  representation  thus  constructed  contains  in  it,  as  parameters,  the  nodal 
values  of  the  function  which  is  being  sought.  To  this  stage  then,  the  problem  is  basically  one  of 
interpolation,  and  the  “interpolation  functions”  over  each  element  are  chosen  so  that  appropriate 
continuity  properties  exist  at  the  boundaries  of  adjacent  elements. 

The  problem  now  has  been  reduced  to  the  determination  of  the  nodal  values  of  the  sought 
function.  This  is  achieved  by  the  application  of  the  relevant  physical  principles  which,  for 
elastic  problems  at  least,  generally  involves  the  use  of  a “minimum  energy"  theorem.  Since  the 
strain  energy  in  an  elastic  body  is  simply  the  sum  of  the  strain  energy  in  its  constituent  parts, 
the  minimisation  of  the  energy  can  be  carried  out  over  the  individual  elements.  Incorporating 
the  results  for  all  elements,  a set  of  simultaneous  equations  is  obtained  for  the  nodal  values  of 
the  function  of  interest.  After  the  application  of  any  boundary  conditions,  these  equations  can 
be  solved  and  the  problem  is  then  essentially  completed. 


1.  Ergatoudis,  I.,  Irons,  B.  M.,  and  Zienkiewicz,  O.  C.,  Curved  Isoparametric  Quadrilateral 
Elements  for  Finite  Element  Analysis.  Inter.  J.  Solids  and  Structures,  vol.  4,  pp.  31-42,  1968. 

2.  Zienkiewicz,  O.  C.,  et  a!..  Isoparametric  and  Associated  Element  Families  foi  Two-  and 
Three-Dimensional  Analysis,  pp.  383-432  of  “Finite  Element  Methods  in  Stress  Analysis", 
edited  by  I.  Holand  and  K.  Bell,  Tapir  Forlag,  Trondheim,  1970. 

3.  Gallagher,  R.  H„  Finite  Element  Analysis  Fundamentals,  Prentice-Hall,  Englewood  Cliffs, 
1975. 


The  above  remarks  apply  to  any  finite  element  problem  and,  in  principle,  elements  of  any 
shape  can  be  used.  However,  the  practical  difficulties  associated  with  the  choice  of  the  interpo- 
lation functions  and  with  the  energy  minimisation,  which  generally  involves  an  integration  over 
the  volume  of  each  element,  are  much  reduced  when  the  elements  have  simple  shapes;  for  example, 
in  two-dimensional  problems  rectangular  or  triangular  elements  are  commonly  used.  The 
difficulty  here  is  that  it  may  sometimes  be  necessary  to  use  a large  number  of  such  elements  in 
order  to  give  an  adequate  geometrical  representation  of  the  actual  region;  this  is  particularly  so 
when  the  region  has  curved  boundaries.  It  is  under  these  circumstances  that  the  use  of  isopara- 
metric elements  can  be  an  advantage. 

The  first,  and  major,  concept  with  regard  to  isoparametric  elements  is  that  an  element  of 
more  or  less  arbitrary  shape  can  be  conveniently  utilised  by  “mapping"  it  on  to  a simpler  shaped 
region  (most  commonly,  a square);  all  subsequent  calculations  are  then  carried  out  for  this 
simple  region. 

The  second  concept  is  that  the  same  functional  form  that  is  used  to  effect  the  mapping  can 
also  be  used  as  the  interpolation  function  for  the  element.  It  is  this  use  of  the  same  function 
for  the  mapping  and  the  interpolation,  which  gives  rise  to  the  name  “isoparametric”.  Although 
there  is  no  fundamental  reason  why  the  same  functional  form  need  be  used  for  both  purposes, 
nevertheless  it  is  often  the  most  convenient  arrangement  in  practice;  see,  for  instance,  p.  410 
of  Reference  2. 


2.2  Interpolation  Functions  for  Square  Elements 

The  following  discussion  relates  to  two-dimensional  problems  which  are  to  be  treated  using 
general  “quadrilateral  elements”.  (The  sides  of  the  elements  may  be  straight  or  curved.)  For 
the  present  purposes,  any  quadrilateral  element  in  the  actual  xy  plane  will  be  mapped  on  to  the 
square  region  defined  by 


— I « ( < 1 


— I ^ v < 1 


} (2.2.1) 


in  the  transformed  fij  plane;  see  Figure  I.  However,  before  discussing  details  of  the  mapping 
it  will  be  useful  to  consider  interpolation  functions  for  the  transformed  region. 


2.2.1  Interpolation  Functions  for  Four-Node  Square 

Suppose  the  square  region  in  the  transformed  plane  has  nodes  at  each  of  its  corners,  identi- 
fied by  the  Roman  numerals  I to  IV  in  Figure  2.  Let  ^<f,  t ;)  be  some  function  which  is  to  be 
expressed  throughout  this  region  in  terms  of  its  values  'Yt  (/  = I to  IV)  at  the  nodes.  It  can  be 
readily  verified  that  this  is  achieved  by  the  expression 

,-?)  = I P<((,  vWt  (2.2.2) 

i - 1 


where 

/Mf,  *?)  = (!  + fXl  +q)/4 
P\\(t  »j)  = (l  — fXl  + i)/4 
Pm((,  v)  = (1  ~ fXl  - ?)/4 
Piv(f,  v)  = (l  + fXl  -q)/4 


(2.2.3) 


At  node  i,  /*<  = 1 whilst  all  the  other  Ps  are  zero.  For  example,  at  node  I,  f — ij  = 1 and  so 
Pi  = 1 whilst  Pii  = Pui  = Piv  = 0;  hence  Equation  (2.2.2)  returns  the  result  </>  — T'i  as 
required.  Actually,  Equation  (2.2.2)  is  simply  the  standard  formula  for  linear  interpolation 
over  a square  region4;  it  implies  a linear  variation  of  the  function  <f>  along  any  straight  line 
parallel  to  a side  of  the  square. 


4.  McCormick,  J.  M.,  and  Salvadori,  M.  G.,  Numerical  Methods  in  Fortran,  p.  1 19,  Prentice- 
Hall,  Englewood  Cliffs,  1964. 


(a)  Physical  plane 


(b)  Transformed  plane 


FIG.  3 EIGHT  NODE  SQUARE 
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At  a later  stage,  the  partial  derivatives  of  if>  with  respect  to  £ and  7?  will  be  required : 

4>i  = {(1  + iK'Fi  - V„)  - (1  - vXViii  ~ Yiv)}/4 

= {(1  + fK^i  - ^iv)  + (1  - OC^ii  - Tra)}/4 

(Note  that  here  and  throughout,  partial  derivatives  with  respect  to  f and  y will  be  written  using 
the  subscript  notation.) 


} (2.2.4) 


2.2.2  Interpolation  Functions  for  Eight-Node  Square 

Now  consider  the  case  where  the  square  in  the  transformed  plane  has  nodes  at  the  mid- 
points of  each  of  its  sides,  as  well  as  at  its  corners;  the  midpoint  nodes  are  numbered  V to  VIII 
as  shown  in  Figure  3.  In  this  case,  the  value  of  </i  at  any  interior  point  can  be  expressed  in  terms 
of  its  nodal  values  VF)  (i  = I to  VIII)  by  the  relation 
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(2.2.5) 


(2.2.6) 


Again,  it  is  readily  established  that,  at  node  i,  Qi  = 1 whilst  all  the  other  Qs  are  zero,  e.g.  at 
node  VI,  f = — 1 and  17  = 0 so  Qvi  = 1,  with  all  other  Qs  zero.  Equation  (2.2.5)  implies  a 
quadratic  variation  of  ip  along  any  straight  line  parallel  to  a side  of  the  square.  However,  it  should 
be  pointed  out  that  Equation  (2.2.5)  is  not  the  standard  quadratic  interpolation  formula  for  a 
square  region,  since  this  last  also  requires  incorporation  of  a nodal  value  at  the  centre  of  a 
square  (Reference  4). 

Analogously  to  Equation  (2.2.4),  the  partial  derivatives  of  Equation  (2.2.5)  are 


<h  = [(I  + vW(  + v)Vi  + (2f  - r,yrn}  + (1  - v){(2£  + O'Fm  + (2f  - #iv} 
- 4f{(l  + i))Tv  + (1  ~ rrt'Tvii}  - 2(1  - VOOFvi  - 'Fvm)]/4 
= [(1  + 0{(f  + 2tj)'Fi  + (-f  + 2r,)>F,v}  + (1  - m-S  + 27?)rII+(f  + 2,)'F,1I} 
+ 2(1  - f*)0Fv  - Tvn)  - 4,{(1  - flTvi  + (1  + fJ'F viii}]/4 


(2.2.7) 


2.3  Mapping  Functions 


Any  pair  of  functional  relations  of  the  form 


X = x(t,  rj) 

y = At  0 


(2.3.1) 


can  be  interpreted  as  a mapping  of  some  region  in  the  x y plane  on  to  some  other  region  in  the 
plane.  The  properties  of  such  mappings  are  discussed  at  length  in  many  mathematical  texts, 
e.g.  that  by  Courant.5  Before  going  on  to  consider  specific  mapping  functions,  some  of  the  general 
results  that  will  be  required  throughout  will  be  set  down. 


5.  Courant,  R„  Differential  and  Integral  Calculus,  vol.  II,  pp.  133  et  seq.,  Blackie,  London, 
1936. 


1 
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If  T(x,  y)  is  a function  defined  over  some  region  R of  the  xy  plane,  then  on  substituting  for 
x and  y from  Equation  (2.3.1),  there  is  obtained  a function  7?)  defined  over  some  other  region 
r of  the  (r)  plane.  The  following  formulae  relate  the  values  of  integrals  and  derivatives  in  the 
actual  plane  to  those  in  the  transformed  plane.  For  integrals, 

SS  Y(x,  y)  dx  dy  = //  V)J((,  v)dfdv  (2.3.2) 

R r 

where  J(£,  17)  is  the  Jacobian  of  the  mapping  and  is  given  by 

J = ~ x,y(  (2.3.3) 

For  derivatives, 

'Vx  = </>(£x  + <\>f)x 

= ifi  + <I>VT)V  (2.3.4) 

In  their  present  form,  the  evaluation  of  the  right  hand  sides  of  Equations  (2.3.4)  requires  a 
knowledge  of  the  inverse  mapping,  i.e.  requires  the  relations 

i = ((x,  y)  (2.3.5) 

v = v(x,  y) 

However,  the  need  for  inverting  Equations  (2.3.1)  to  obtain  Equations  (2.3.5) — which  is  often 
not  convenient  in  practice — can  be  obviated  by  the  use  of  the  relations 

£x  = yJJ  = -xJJ  (2.3.6) 

Vx  — —yJJ  Vu  — XfjJ 

where  J is,  as  before,  the  Jacobian.  Hence,  substituting  from  (2.3.6)  into  (2.3.4)  gives 

^ = (0£y,  -MIJ  (2.3.7) 

v*  = (~hxv  + MU 

The  evaluation  of  the  right  hand  sides  of  Equations  (2.3.7)  only  requires  a knowledge  of  the 
original  mapping  (2.3.1). 

Now  attention  will  be  turned  to  specific  mapping  functions  that  are  utilised  in  the  case  of 
isoparametric  quadrilateral  elements. 

2.3.1  Mapping  of  Straight-Sided  Quadrilateral  on  to  Square 

Consider  an  arbitrary  straight-sided  quadrilateral  in  the  xy  plane  and  let  xi,  yi  (i  = I to  IV) 
denote  the  co-ordinates  of  the  comers  (Fig.  4).  Define  a mapping  by 


x = £ Pt(L  v)xt 

t = 1 

y = £ Pt((,  v)yi 

i = 1 


(2.3.8) 


where  the  Pi  are  given  by  Equations  (2.2.3).  This  mapping  transforms  exactly  the  quadrilateral 
into  the  square  region  of  Figure  2 in  the  transformed  plane.  This  can  be  seen  as  follows.  Each 
corner  of  the  quadrilateral  certainly  maps  into  the  corresponding  corner  of  the  square  since 
Pi  = 1 at  the  corner  “i”  of  the  square  and  all  the  other  P s are  zero,  so  Equations  (2.3.8)  simply 
return  the  result  x = xi,  y = yt.  Further,  a side  of  the  quadrilateral  maps  exactly  into  a side  of 
the  square.  For  example,  the  line  joining  corners  I and  II  of  the  square  corresponds  to  17  = 1 
in  the  transformed  plane.  On  substituting  this  value  into  Equations  (2.3.8),  these  last  reduce  to 


x = {(1  + f)xi  + (1  — f)*n}/2 

y = {(I  + ()yi  + ( 1 - ()yn}/2 

On  eliminating  f between  these  equations  the  result  is 

Ax  11  — *1)  = x(yn  - yt)  + xuyi  — xtyu 


(2.3.9) 


(2.3.10) 


NODAL  COORDINATES  FOR  STRAIGHT-SIDED  QUADRILATERAL 


I 


NODAL  COORDINATES  FOR  CURVILINEAR  QUADRILATERAL 


This  is  the  equation  of  the  straight  line  joining  the  corners  1 and  II  of  the  quadrilateral.  Similar 
results  hold  for  all  sides  and,  hence,  Equations  (2.3.8)  map  the  sides  of  the  quadrilateral  exactly 
into  the  sides  of  the  square. 

The  various  derivatives  needed  for  evaluating  the  Jacobian  (2.3.3)  and  which  are  also  needed 
in  Equations  (2.3.7)  can  be  readily  obtained  from  Equation  (2.3.8): 


x(  = {(1  + t/Xxi  — xu)  — (1  — tjXxiii  — *iv)}/4 

x,  = {(1  -f  f)(xi  — xiv)  + (1  — £)(*n  — Xni)}/4 

>’£  = {(•  + v)(yi  ~ ^ii)  — (1  — ’/XJ'iii  — yiv)}/4 

yv  = {(1  + f)(>’i  — >'iv)  + (1  — fXj'n  — yni)}/4 


(2.3.11) 


2.3.2  Mapping  of  Curvilinear  Quadrilateral  on  to  Square 

Now  consider  a curvilinear  quadrilateral  in  the  xy  plane  and  let  xi,  yi  (i  = I to  IV)  denote 
the  co-ordinates  of  the  corners  and  xt,  yt  (/  = V to  VIII)  denote  the  co-ordinates  of  intermediate 
points  on  the  sides  of  the  quadrilateral,  ordered  as  shown  in  Figure  5.  No  stipulation  is  made 
about  particular  locations  for  the  intermediate  nodes  V to  VIII.  Define  a mapping  by 


VIII 

x = E Qt(L  v)xt 

t = i 

VIII 

y = L (?<(£  v)yt 


(2.3.12) 


where  the  Qs  are  given  by  Equations  (2.2.6).  This  mapping  transforms  approximately  the  curvi- 
linear quadrilateral  in  the  xy  plane  into  the  square  of  Figure  3 in  the  £77  plane.  By  virtue  of  the 
properties  of  the  Qs,  Equations  (2.3.12)  will  map  the  eight  “nodal”  points  in  the  xy  plane 
exactly  into  the  corresponding  eight  nodal  points  in  the  £77  plane.  However,  in  distinction  to  the 
situation  for  a straight-sided  quadrilateral,  for  non-nodal  points  on  the  boundary  of  the  square 
in  the  fy  plane,  the  values  of  x and  y as  computed  from  Equations  (2.3.12)  will  not  in  general  lie 
on  the  boundary  of  the  curvilinear  quadrilateral,  i.e.  the  mapping  is  not  an  exact  one.  (Of  course, 
this  could  not  be  expected  since  the  quadrilateral  has  only  been  specified  to  the  extent  of  its  eight 
nodes.  Whilst  these  nodes  have  uniquely  specified  the  mapping  function  (2.3.12),  they  certainly 
have  not  uniquely  specified  the  quadrilateral.) 

The  general  nature  of  the  quadrilateral  that  is  being  mapped  by  (2.3.12)  can  be  established 
by  examining  what  happens  as  one  traverses  a side  of  the  square  in  the  £77  plane.  For  example, 
referring  to  Figure  3 it  can  be  seen  that  the  side  I-V-II  corresponds  to  77  = 1 and,  on  substituting 
this  into  (2.3.12)  the  result  is 


x = HI  + £)*i/2  - |(1  - 0*ii/2  + (1  - 0)*v 
y = f(l  + 0W2  - 01  - Ojdi/2  + (I  - P)yv 

When  f is  eliminated  from  these  two  equations,  a result  of  the  form 


(2.3.13) 


fix2  + wy  + C3>’2  + c\x  + csy  + d = 0 (2.3.14) 

is  obtained,  where  the  coefficients  ry  are  rather  cumbersome  functions  of  the  co-ordinates 
xt,yt  (/ = I,  II,  V).  Equation  (2.3.14),  which  represents  a second  degree  curve,  defines  the 
boundary  curve  in  the  xy  plane  which  is  actually  being  mapped.  So  far,  node  V has  not  been 
specified  in  the  xy  plane  beyond  the  requirement  that  it  lie  somewhere  along  the  side  I— II.  If 
the  choice 

*v  = (*i  + *n)/2  (2.3.15) 

be  made,  then  the  first  of  Equations  (2.3.13)  simplifies  to 

x = {((xi  — xn)  + (xi  + xn)}/2  (2.3.16) 

On  eliminating  f between  (2.3.16)  and  the  second  of  (2.3.13)  the  result  is  of  the  form 

y = </ix2  + diX  + ds  (2.3.17) 
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where  the  coefficients  dj  are  functions  of  the  .v<  and  yi.  Here,  then,  the  boundary  curve  in  the  x y 
plane  is  simply  the  (unique)  parabola  through  the  three  nodal  points.  In  this  case  the  accuracy 
of  the  mapping  is  governed  by  the  closeness  with  which  the  actual  boundary  curve  can  be 
approximated  by  a parabola. 

However,  it  is  not  necessary  to  make  the  choice  (2.3.15)  and,  in  any  given  case,  the  map- 
ping actually  being  achieved  by  Equations  (2.3.12)  can  be  established  by  a direct  evaluation.  As 
an  example,  consider  the  quadrant  of  a circle  of  unit  radius  shown  in  Figure  6.  Later  on,  the 
torsion  problem  for  a bar  of  this  cross-section  will  be  worked  out  in  detail  using  the  (coarse) 
finite  element  subdivision  shown.  For  the  moment  attention  is  restricted  to  element  “b”.  The 
nodal  numbering  along  with  the  nodal  co-ordinates  for  this  element  are  shown  in  Figure  7. 
The  co-ordinates  of  the  corner  nodes  are,  of  course,  determined  by  the  finite  element  subdivision. 
The  intermediate  nodes  have  been  selected  so  as  to  occur  half  way  along  the  sides.  (In  particular, 
since  the  arc  IV-I  subtends  an  angle  of  u-/4  at  the  origin,  the  arc  IV— VIII  subtends  an  angle  of 
w/8;  the  relation  analogous  to  (2.3.15)  does  not  apply  here.)  On  substituting  the  values  of  x<  and 
>’(  as  given  on  Figure  7 into  Equations  (2.3.12)  the  mapping  function  for  element  “b"  is  defined. 
In  order  to  assess  the  accuracy  of  the  mapping,  values  of  x and  y will  be  calculated  from  (2.3.12) 
for  various  values  of  f and  i)  corresponding  to  points  on  the  boundary  of  the  square  of  figure  3. 
Attention  can  be  restricted  to  the  curved  side  IV— VIII-I  because  the  mapping  is  exact  for  any 
straight  side,  as  can  easily  be  proven.  The  side  IV— VIII— I corresponds  to 

f = 1,  -1  < i?  « 1 (2.3.18) 

in  the  transformed  plane.  On  substituting  f = 1 into  (2.3.12)  and  inserting  there  the  values  of 
xt  and  yt  as  obtained  from  Figure  7,  the  mapping  function  reduces  to 

x = 0-7071i7(l  + 7?)/2  - 1 -000017(1  - i;)/2  + 0-9239(1  - i,2)  (2.3.19) 

y = 0-7071i7(l  + rj)/2  + 0 + 0-3827(1  - i)2) 

The  result  of  evaluating  (2.3.19)  for  various  values  of  ij  is  shown  in  Table  1 below.  As  a measure 
of  the  error  in  the  mapping,  the  difference,  e,  between  the  radius  vector  from  the  origin  to  the 
mapped  point  and  the  radius  of  the  circular  boundary  (unity)  has  also  been  tabulated.  Here, 

c = (x2  + y2)*  - 1 (2.3.20) 

It  can  be  seen  that  the  error  is  everywhere  small  despite  the  relatively  large  length  of  arc 
involved.  (The  error,  of  course,  is  zero  at  the  nodal  points.) 


| 


TABLE  1 

Evaluation  of  Mapping  Function  for  Circular  Arc 


V 

X 

y 

€ 

-1-00 

1 -0000 

0-0000 

0-0000 

-0-75 

0-9942 

0-1011 

-0-0007 

-0-50 

0-9795 

0-1986 

-0-0005 

-0-25 

0-9561 

0-2925 

-0-0002 

0-00 

0-9239 

0-3827 

0-0000 

0-25 

0-8829 

0-4693 

-0-0002 

0-50 

0-8331 

0-5522 

-0-0005 

0-75 

0-7745 

0-6315 

-0-0007 

1 -00 

0-7071 

0-7071 

0-0000 

Finally,  the  various  derivatives  needed  for  use  in  Equations  (2.3.3)  and  (2.3.7)  will  be  set 
down: 
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xf  = [(1  + + y)xi  + (2f  — ij)jfn}  + (1  — r;){(2f  + >y)jrni  + (2{  — i))xiv} 

— 4£{(  l + rj)xv  + ( 1 — i?)xvn}  —2(1—  i?2X-*'vi  — Jfvm)]/4 

(2.3.21) 

x,  = [(1  + fXtf  + 2^  + (-(  + 2v)x,v}  + (1  - (X(-(  + 2r,)xu  + ((  + 2r,)xm} 

+ 2(1  - PHxv  ~ xvn)  - 4tj{(  1 - f)xvi  + (1  + 0*vm}]/4 
The  formulae  for  j>{  and  are  obtained  simply  by  replacing  x,  by  y,  (/'  = 1 to  VIII)  in  the  above. 

3.  FORMULATION  OF  TORSION  PROBLEM  USING  ISOPARAMETRIC 
ELEMENTS 

3.1  General 

In  the  next  two  Sections  the  use  of  isoparametric  elements  will  be  demonstrated  by  examples 
involving  the  torsion  of  bars  and  it  is  convenient  to  set  down  here  the  basic  equations  governing 
the  torsion  problem.  Consider,  therefore,  a bar  of  (constant)  cross-section,  R,  and  let  c denote 
the  boundary  curve  of  the  cross-section.  It  is  shown  in  standard  texts  on  elasticity,  e.g.  that  by 
Sokolnikoff,8  that  the  solution  of  the  torsion  problem  can  be  reduced  to  finding  a function 
TX*,  y)  such  that 


'Vxx  + T'VJ,  = —2  over  R 
T"  = 0 on  c 


(3.1.1) 

(3.1.2) 


where  x and  y are  the  co-ordinates  in  the  plane  of  the  cross-section  and  the  suffixes  again  denote 
partial  differentiation.  Once  'F  is  determined  the  torsion  constant,  K,  is  given  by 

K=2jS'V(x,y)dxdy  (3.1.3) 

R 

and  the  shear  stresses  rx  and  ry  by 

TzlT='YylK , tv/T=-'Yx/K  (3.1.4) 

where  T denotes  the  applied  torque.  The  angle  of  twist  per  unit  length,  6,  is  related  to  the  torque 
according  to 


r=  GKO 


(3.1.5) 


where  G is  the  shear  modulus. 

3.2  Integral  Formulation 

For  finite  element  work,  the  differential  equation  formulation  embodied  in  Equation  (3.1.1) 
is  not  suitable.  Instead,  use  is  made  of  the  fact  that  the  function  V which  satisfies  (3.1.1)  also 
makes  the  integral 

/(Y)  = (1/2)  JJ  (MV  + TV  - 4 V)dxdy  (3.2.1) 

R 

a minimum.  This  equivalence  can  be  established  either  purely  mathematically  using  the 
Calculus  of  Variations  or,  physically,  by  using  the  principle  of  Minimum  - Jomplementary  Energy 
(see  p.  416  of  Ref.  6).  The  boundary  condition  (3.1.2),  of  course,  still  rr.  n be  satisfied. 

The  minimisation  required  in  (3.2.1)  is  generally  carried  out  as  follows.  It  is  assumed  that 
Y can  be  written  in  the  form 


%;)  = ! fk{x,y) Y* 

k = I 


(3.2.2) 


where  the  /*  are  known  functions,  and  are  such  that  the  boundary  condition  (3.1.2)  is  satisfied, 
whilst  the  Y*  are  presently  unknown  parameters.  (As  will  soon  transpire,  the  Y*  will  be  identi- 

6.  Sokolnikoff,  I.  S„  Mathematical  Theory  of  Elasticity,  2nd  edition,  p.  116,  McGraw-Hill, 
New  York,  1956. 


10 


fied  with  the  nodal  values  of  Y in  a finite  element  subdivision  of  the  region  R.)  On  substituting 
from  (3.2.2)  into  (3.2.1)  it  can  be  seen  that  the  integral  I becomes  a function  of  the  Y*,  i.e. 


/(Y)  = /(Y,, . . . Y,,) 

The  minimisation  of  / is  achieved  by  solving  the  n simultaneous  equations 


*-.0 

JY* 


k = 1,. . .,/» 


(3.2.3) 


(3.2.4) 


for  the  Y*.  Once  these  have  been  determined  the  problem  is  basically  solved.  (Note  that  partial 
differentiation  with  respect  to  the  Y*  will  always  be  denoted  by  the  “3”  symbol,  rather  than  a 
suffix;  however,  the  suffix  notation  will  continue  to  be  used  for  partial  differentiation  with 
respect  to  x,  y,  f and  ij.)  The  above  procedure  is  basically  an  approximate  one,  but,  in  general, 
the  approximation  can  be  rendered  adequate  by  taking  a sufficient  number  of  Y*. 

3.3  Finite  Element  Formulation 

In  the  finite  element  approach  to  the  torsion  of  a bar,  the  cross-section  R is  considered  as 
being  subdivided  into  a number  of  elements  Ri.  (A  simple  example  has  already  been  shown 
in  Fig.  6.)  Because  of  the  additive  nature  of  integrals  it  is  possible  to  write,  in  place  of  Equation 
(3.2.1) 

1=1.1, 1 (3-3.1) 

/ 


U = (1/2)JJ(Y,2  + Y„2  - 4Y )dxdy 


Then  Equations  (3.2.4)  become 


1 — = 0 

, 3Y* 


k = 


(3.3.2) 


(3.3.3) 


Hence,  attention  is  concentrated  initially  on  evaluating  the  derivatives  3/j/iY*  for  each  element; 
then  the  final  set  of  equations  is  obtained  by  carrying  out  the  summation  (3.3.3).  From  Equation 
(3.3.2) 


3 = IS  U v W,  _ 2^-\dxdy 

3Y*  *}  vSy*  + 4|,3Y*  23Y kf  y 


(3.3.4) 


In  a direct  application  of  the  finite  element  method,  some  form  of  variation  of  Y over  each 
element  is  now  assumed,  and  the  integral  (3.3.4)  is  evaluated.  However,  as  already  mentioned, 
there  are  practical  difficulties  for  any  but  the  simplest  shaped  elements,  firstly,  in  choosing  an 
appropriate  form  for  Y and,  secondly,  in  actually  evaluating  the  integral. 

Before  passing  on  it  might  be  noted  that  the  torsion  constant  K can  be  written  in  the  form 

K = 'ZKi  (3-3.5) 


Ki  = 2 Jf  Y (x,y)dxdy 

» i 


(3.3.6) 


3.4  Isoparametric  Element  Formulation 

Suppose  that  the  subdivision  of  the  cross-section  R is  made  using  either  the  four-node 
elements  of  Figure  4 or  the  eight-node  elements  of  Figure  5.  Then,  in  order  to  obviate  the 
difficulties  just  mentioned,  each  element  Ri  in  the  xy  plane  is,  in  turn,  mapped  on  to  the  square 
region  in  the  fr;  plane  either  by  Equation  (2.3.8)  or  by  Equation  (2.3.12)  as  appropriate.  In  this 
case  the  function  Y(jr,  y)  will  transform  to  a function  </<f,  tj)  and,  using  the  relations  (2.3.2) 
and  (2.3.7),  Equation  (3.3.4)  becomes 


a/i 

i'Vi 


+ (-*■«. + + ?t,*‘)\/j 
- «£]«* 


(3.4.1) 


Note  that,  since  attention  is  at  present  being  confined  to  a single  element  (the  “/th”),  a local 
nodal  numbering  system  can  be  used  temporarily.  Hence,  T<  has  been  written  in  (3.4.1)  in  place 
of  T*  in  (3.3.4)  where  i takes  the  values  I to  IV  for  a four-node  element  and  I to  VIII  for  an 
eight-node  element.  (Of  course,  at  a later  stage  the  local  numbering  must  be  replaced  by  a global 
numbering.) 

Equation  (3.4.1)  can  be  written  in  the  form 


where 


A = x*  + yf 

B = xtxi  + yty*  (3.4.3) 

C = V + yt2 


Whilst  Equation  (3.4.2)  may  appear  more  formidable  than  Equation  (3.3.4)  its  evaluation  is 
straightforward,  albeit  tedious.  Depending  on  the  type  of  element  being  used,  4>  is  taken  either 
in  the  form  (2.2.2)  or  (2.2.5),  with  <fi(  and  </>,  correspondingly  being  given  by  either  (2.2.4)  or 
(2.2.7).  In  both  cases  one  can  write 


j 

'/',  = ! v,r,  (3.4.4) 

J 

4>  = Y 

j 

where  the  us,  vs  and  ws  are  functions  of  £ and  r)  whose  explicit  forms  can  be  obtained  from  the 
equations  just  cited;  these  are  listed  in  Tables  2 and  3 below.  Clearly, 

u(  = 2^/JVi,  v,  = W,IWu  wi  = ty/S'F<  (3.4.5) 

In  Equations  (3.4.4)  the  summation  goes  from  j = I to  IV  for  a four-node  element  and  from 
j = I to  VIII  for  an  eight-node  element.  On  substituting  from  (3.4.4)  into  (3.4.2)  it  follows  that 

^ = (3.4.6) 

J 

where 

a<j  = J jJ  , ((^Ui  ~ Bv<)ul  + (Cvt  ~ But)vf}/Jd(dT)  (3.4.7) 

*<  = 2 J1  J1 1 WiJdidr,  (3.4.8) 

The  integrations  (3.4.7)  and  (3.4.8)  are  virtually  always  done  numerically  using  a Gaussian 
integration  formula.7 


7.  Hildebrand,  F.  B.,  Introduction  to  Numerical  Analysis,  McGraw-Hill,  New  York,  1956. 
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TABLE  2 

Formulae  for  u<,  v(,  w<  for  Four-Node  Element 


1 

Ut 

V| 

I 

(1  + i))/4 

(1  + 0/4 

(1  + 0(1  + ij)/4 

II 

-(i  + V)l* 

(i  - 0/4 

(1  - 0(1  + ij)/4 

III 

-(1  - ij)/4 

-0  - 0/4 

(1  — 00  — i?)/4 

IV 

(1  - 1»)/4 

-0  + 0/4 

(1  + 00  -*)/4 

TABLE  3 

Formulae  for  u<,  v<,  w<  for  Eight-Node  Element 


Ui 

0 

+ VX2(  + rf)/4 

(1 

+ .»X2f  - 1?)/4 

(1 

- ,X2f  + i?)/4 

0 

- VX2(  - ff)/4 

-«i  + v) 

-0  - v*)l2 

-«i  -V) 

(1  - r/2)/2 

(1  + fXf  + 2t,)/4 
(1  -(X-t  + 2v)l4 
(i  - m + 2i?)/4 
(1  + fK-f  + 2tj)/4 
(1  - **)/2 
-iKl  - 0 
-(1  - f*)/2 


(1  + fl(l  + ,x-l  + ( + v)/4 
(1  -0(1  +1JX-1  — f + i»)/4 
(1  -00  -iX-1  -i-v)l4 
(1  + 0(1  - i»X-l  + ( ~ v)l* 
(i  + ,xi  - a/2 

(1  _ fKi  _ ,« )/2 
(i  _,xi  - a/2 

a + 0(1  - a/2 


When  the  quantities  and  have  been  determined  for  each  element,  the  global  equations 
(3.3.3)  can  be  set  up.  The  boundary  condition  (3.1.2)  is  then  applied  by  setting  T*  = 0 for 
each  boundary  node.  After  solving  the  resulting  set  of  equations  for  the  remaining  T*  the 
problem  is  essentially  complete.  The  torsion  constant  K can  be  computed  from  Equation  (3.3.5) 
where  now,  in  place  of  Equation  (3.3.6) 


<K(,  *?)•/( f.  v)di(h) 


(3.4.9) 


On  substituting  from  the  third  of  Equations  (3.4.4)  into  Equation  (3.4.9),  and  using  Equation 
(3.4.8),  it  follows  that 

= (3.4.10) 

the  summation  going  either  from  I to  IV  or  I to  VIII. 

The  stresses  can  be  obtained  by  substituting  from  Equations  (2.3.7)  into  Equations  (3.1.4) 
to  get 

rz/T  = (-*,*,  + M/(JK)  (3.4. 1 1) 

r,/T= 

4.  ILLUSTRATIVE  EXAMPLE  USING  FOUR-NODE  ELEMENTS— TORSION  OF 
BAR  WHOSE  SECTION  IS  A RIGHT  ANGLE  ISOSCELES  TRIANGLE 

4.1  General 

As  an  example  of  the  use  of  four-node  isoparametric  elements,  the  torsion  of  a bar  having  a 
cross-section  in  the  form  of  a right-angle  isosceles  triangle  will  be  considered.  The  vertices  of 
the  triangle  are  located  at  the  points  (1,0),  (0,  1),  (—1,0)  as  shown  in  Figure  8.  The  (crude) 
finite  element  subdivision  adopted  is  shown  in  Figure  9,  the  elements  being  identified  as  “a”, 


n 

i I 


“b”  and  “c".  Also  shown  in  Figure  9 are  the  global  node  numbers  (I  to  7)  and  the  local  node 
numbers  for  each  element  (I  to  IV).  These  last  must  run  sequentially  in  a counter-clockwise 
sense  within  any  one  element;  however,  node  I can  be  assigned  arbitrarily.  The  co-ordinates  of 
the  nodes  are  as  shown  in  Table  4 below. 

TABLE  4 

Nodal  Co-ordinates  for  Triangle 


Element 

node 

no. 


Element  a 

Element  b 

Element  c 

Global 

Global 

Global 

node 

Xi 

yt 

node 

Xi 

yt 

node 

yt 

no. 

no. 

no. 

4 

0 

0-3333 

4 

0 

0-3333 

4 

0 

0-3333 

1 

o 

0 

5 

-0-5 

0-5 

3 

0-5 

0-5 

2 

1 

0 

6 

-1 

0 

7 

0 

1 

3 

0-5 

0-5 

1 

0 

0 

5 

-0-5 

0-5 

Table  4 contains  all  the  data  needed  for  the  construction  of  the  mapping  functions  and  their 
derivatives.  However,  before  passing  on  to  these,  some  reference  should  be  made  to  the  numerical 
integration  scheme  which  will  be  used  throughout  this  Section  and  the  next.  The  determination  of 
the  at)  and  which  are  defined  by  Equations  (3.4.7)  and  (3.4.8)  always  requires  the  evaluation 
of  an  integral  of  the  form 


(4.1.1) 


Here,  the  four-point  Gaussian  formula,  namely, 

F*A~  1/V3,  -l/\/3)  +/(I/a/3,  -1/V3)  +/(-l/V3,  l/y/3)  +/(I/v/3,  1/^3) 

(4.1.2) 

will  always  be  used.  In  practice,  a higher  order  formula,  e.g.  the  nine-point  formula,  should 
generally  be  used  but  (4.1.2)  will  serve  to  illustrate  the  nature  of  the  calculations. 

4.2  Mapping  Functions 

Actually,  the  mapping  functions  for  the  various  elements  are  not  required  anywhere  in 
the  calculations;  it  is  only  the  derivatives  of  the  mapping  functions  which  appear.  However, 
for  completeness  they  will  be  included  here.  Each  element  will  be  considered  in  turn. 

Element  a 

Substituting  the  values  of  x<  and  y<  (/  = I to  IV)  from  Table  4 into  Equations  (2.3.8)  gives 


* = {(1  + #xi  + v)(0)  + (l  - 0(1  + 0(0)  + (i  - 0(1  - 0(1) 


+ (1+0(1  -O(0-5)}/4 

y = {(1  + 0(1  + 0(0-3333)  + (1  - 0(1  + 0(0)  + (1  - 0(1  - ijXO) 

+ (1  +0(1  -O(0-5)}/4 

These  simplify  to 

x = 0-375  -0-125f  - 0-375i,  + 0-1250)  (4.2.1) 

y = 0-2083  + 0-2083f  - 0-0416i,  - 0-04160, 

The  derivatives  of  the  mapping  function  as  calculated  from  Equations  (2.3.11)  (or  from  (4.2.1)) 


x(  = -0-125(1  -0  x,  = -0  125(3  - 0 
^ = 0-0416(5-1,)  y,  = -0-0416(1  + 0 


(4.2.2) 


(It  might  be  observed  that  all  calculations  were  done  using  more  figures  than  will  be  displayed 
here;  the  displayed  numbers  here  and  throughout  have  been  obtained  by  truncating,  rather  than 
rounding,  the  actual  numbers.  This  gives  rise  to  some  apparent  discrepancies  of  a minor  nature.) 

It  is  now  necessary  to  evaluate  the  derivatives  (4.2.2)  for  values  of  ( and  i,  corresponding  to 
the  points  utilised  in  the  Gaussian  integration.  The  results  are  shown  in  Table  5. 


TABLE  5 

Derivatives  of  Mapping  Function  for  Element  a 


The  values  of  Table  5 are,  in  turn,  used  to  calculate  the  quantities  J,  from  Equation  (2.3.3), 
and  A,  B,  C from  Equations  (3.4.3).  The  results  are  shown  in  Table  6. 


TABLE  6 

Values  of  J,  A,  B,  C for  Element  a 


V 

-1V3,  — I/V3 

1/V3,  -1/V3 

-l/v'3,  1/V3 

1/V3,  l/v/3 

J 

01073 

0-0833 

0-0833 

0-0592 

A 

0-2002 

0 0960 

0-2002 

0-0960 

B 

0 0840 

0-0445 

0-0203 

0-0038 

C 

0-0928 

0-0928 

0-0367 

0-0367 

Element  b 

There  is  no  need  to  carry  out  any  calculations  for  element  b,  because  the  final  result  for  it 
can  be  determined  from  that  for  element  a by  appealing  to  symmetry. 

Element  c 

Again  the  mapping  function  will  be  displayed.  Using  the  values  from  Table  4,  appropriate 
to  element  c,  this  is  found  to  be 

x = 0-25(-f  + ,)  (4.2.3) 

y = 0 0833(7  - 2(  - 2V  + (v) 

The  derivatives  of  the  mapping  function  as  calculated  from  Equations  (2.3.11)  (or  (4.2.3))  are 

x(  — -0-25  x,  = 0-25  (4.2.4) 

y<  = 0-0833(-2  + i,)  y,  = 0-0833(-2  + () 

The  values  of  these  derivatives  at  the  points  used  in  the  Gaussian  integration  are  shown  in 
Table  7 and  the  subsequently  calculated  values  of  the  quantities  J,  A,  B and  C are  shown  in 
Table  8. 
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TABLE  7 

Derivatives  of  Mapping  Function  for  Element  c 


(>v 

-l/v/3.  -l/v/3 

1/V3,  -l/v/3 

-l/v/3,  l/v/3 

1/V3,  l/v/3 

-0-2500 

-0-2500 

-0-2500 

rm 

0-2500 

1 

0-2500 

0-2500 

-0-2147 

Ell 

-0-1185 

H 

-0-2147 

EH 

Bi 

-01 185 

TABLE  8 

Values  of  J,  A,  B,  C for  Element  c 


€,V 

-1/V3,  -1/V3 

1/V3  , -l/v/3 

—l/v/3,  l/v/3 

l/v/3,  l/v/3 

rm 

0-1073 

■ 

MEM 

0-1086 

0-0765 

Bca 

■si 

B 

-0-0163 

-0-0370 

C 

0-1086 

0-1086 

Boa 

■SI 

4.3  Element  Matrices 


The  next  stage  of  the  calculation  is  the  determination  of  the  quantities  an  and  bi,  as  defined 
by  Equations  (3.4.7)  and  (3.4.8)  for  each  element.  The  ay,  of  course,  are  the  components  of  a 
symmetric  4x4  matrix  for  each  element.  However,  it  is  convenient  to  first  evaluate  the 
quantities  w<,  Vi,  and  h<  (defined  in  Table  2)  at  the  points  used  in  the  Gaussian  integration,  since 
these  have  the  same  values  for  all  elements.  The  results  are  shown  in  Tables  9,  10  and  1 1 below. 


i 


TABLE  9 

Value  of  u<  at  Integration  Points 


i 

—l/v/3,  -l/v/3 

l/v/3,  -l/v/3 

— l/v/3,  l/v/3 

l/v/3,  l/v/3 

i 

0-1056 

0-1056 

ii 

-0-1056 

-0-1056 

hi 

-0-3943 

-0-3943 

BffEi 

■EXgfl 

IV 

0-3943 

0-3943 

■DBS 

■SSI 

TABLE  10 

Values  of  v(  at  Integration  Points 
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i 

— l/v/3,  —l/v/3 

l/v/3,  —l/v/3 

-l/v/3,  l/v/3 

l/v/3,  l/v/3 

i 

0-1056 

fna 

ira 

ii 

1 

HfSl 

KO 

hi 

-0-1056 

^E£2I 

IV 

-0-1056 

-0-3943 

-0-3943 

TABLE  11 

Values  of  w(  at  Integration  Points 


n 

-1/V3.  — 1/V3 

1/V3,  -1/V3 

— 1/3/3.  1/V3 

1/V3,  I/>/3 

i 

0 0446 

0-1666 

0-1666 

0-6220 

ii 

01666 

0-0446 

0-6220 

0-1666 

hi 

0-6220 

K2C9 

0-1666 

0-0446 

IV 

0-1666 

0-0446 

0-1666 

For  the  remainder  of  the  calculations  it  is  necessary  to  consider  each  element  individually. 


Element  a 

As  an  intermediate  step,  the  quantities  (Aui  — Bvt ) and  (Cv<  — But)  are  calculated  at  the 
integration  points  using  the  values  given  in  Tables  6,  9 and  10.  The  results  are  shown  in  Tables 
12  and  13  below. 


TABLE  12 

Values  of  (Au<-Bv<)  at  Integration  Points 


tv 

i 

-W3.  -1/V3 

1/V3,  -1/3/3 

—1/3/3.  1/3/3 

1/V3,  1/3/3 

i 

0-0122 

EBI 

0-0768 

H9I 

ii 

-0-0543 

■ ' 

in 

-0  0458 

-0-0097 

IV 

0-0878 

BEH 

0-0116 

TABLE  13 

Values  of  (Cv(-Bu<)  at  Integration  Points 


tv 

i 

-1/3/3.  — 1/3/3 

1/3/3,  -1/3/3 

-1/3/3.  1/3/3 

1/3/3.  1/3/3 

1 

0 0009 

0-0319 

-0-0041 

0-0129 

11 

0-0145 

0-0054 

hi 

0-0077 

-0-0034 

IV 

-0  0429 

-0-0542 

-0  0060 

-00149 
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Recalling  Equation  (4.1.2),  the  a</  and  bt  for  element  a can  now  be  evaluated  using 
Equations  (3.4.7)  and  (3.4.8)  in  conjunction  with  the  numerical  values  given  in  Tables  6,  9,  10, 
11,  12  and  13.  Three  typical  calculations  are  reproduced  below. 

an  =(0  0122  x 0 1056  + 0 0009  x 0- 1056)/(0- 1073) 

+ (-0  0074  x 0 1056  + 0 0319  x 0 • 3943)/(0 • 0833) 

+ (0  0768  x 0-3943  - 0-0041  x 0- 1056)/(0  0833) 

+ (0-0363  x 0-3943  + 0-0129  x 0-3943)/(0-0592) 

= 0-8407 

034  = (-0-0458  X 0-3943  + 0-0034  X 0- 1056)/(0- 1073) 

+ (—0-0331  X 0-3943  - 0-0077  X 0-3943)/(0  0833) 

+ (—0-0131  X 0 1056  + 0-0123  X 0- 1056)/(0 - 0833) 

+ (-0-0097  X 0 1056  + 0-0034  X 0-3943)/(0  0592) 

= -0-3537 

b2  = 2(0-1666  x 0-1073  + 0 0446  x 0-0833  + 0 6220  X 0-0833 
+ 0-1666  x 0-0592)  = 0-1666 

By  proceeding  in  this  way  the  full  set  of  equations  for  element  a is  found  to  be 


S/a/aV i ' 

a/a/aTn 

a/o/aTni 
a/a/aT  iv 


0-8407  -0-5605  -0-2194  -0-0607' 
-0-5605  1-0404  0-1464  -0-6263 

-0-2194  0-1464  0-4266  -0-3537 

-0-0607  -0-6263  -0-3537  1-0404 


0- 1388 
0-1666 
0-1944 
0-1666 


(4.3.1) 


Using  the  correspondence  between  local  nodes  and  global  nodes  given  in  Table  4 these  last 
equations  can  be  written  as 


a/a/aTY 

a/a/aTi 

a/o/ar2 

a/a/aTg 


Element  b 


0-8407  -0-5605  -0-2194  -0-0607 
-0-5605  1-0404  0 1464  -0-6263 

-0-2194  0-1464  0-4266  -0-3537 

-0-0607  -0-6263  -0-3537  1-0404 


0- 1388 
0-1666 
0-1944 
0-1666 


(4.3.2) 


As  already  remarked  the  results  for  element  b can  be  written  down  by  utilising  symmetry. 
Referring  to  Figure  9 it  can  be  seen  that  nodes  4,  1,  2 and  3 in  element  a correspond  respectively 
to  nodes  4,  1,6  and  5 in  element  b.  Hence,  on  making  the  necessary  changes  in  Equation  (4.3.2), 
the  results  for  element  b become 


a/b/aTY 

a/a/aTr 

a/6/ar« 

a/j./a'Fs 


Element  c 


0-8407  -0-5605  —0-2194  -0-0607 


-0-5605  1-0404 

-0-2194  0-1464 


0-1464  -0-6263 
0-4266  -0-3537 


-0-0607  -0-6263  -0-3537  1-0404 


0- 1388 
0-1666 
0-1944 
0-1666 


(4.3.3) 


The  calculations  for  element  c are  quite  analogous  to  those  for  element  a.  The  quantities 
(Aut  — Bv()  and  (Cv<  — But)  need  to  be  recalculated  using  the  values  given  in  Tables  8,  9 and 
10.  The  mi  and  bi  are  then  calculated  as  before.  Only  the  final  result  will  be  cited  here. 


JJ 


(4.3.4) 


“<>/*/a¥4“ 

' 11818  -0-2878  -0-6060  -0-2878“ 

f¥4-| 

“0-1388" 

s/c/a¥3 

-0-2878  0-5252  -0-0959  -0-1414 

¥3 

0-1666 

a/c/s¥7 

-0-6060  -0-0959  0-7979  -0  0959 

¥7 

0-1944 

_a/c/d¥s_ 

_— 0-2878  -0-1414  -0  0959  0-5252_ 

_¥,_ 

_0- I666_ 

4.4  Global  Equations  and  Solution 

It  is  now  possible  to  set  up  the  global  Equations  (3.3.3)  simply  by  the  appropriate  addition 
of  Equations  (4.3.2),  (4.3.3)  and  (4.3.4).  Each  global  equation  is  of  the  form 


5/  2>/o  3/&  3/c 

i¥*  = 5Y*  + ST*  + 3¥* 


(4.4.1) 


For  example,  the  equation  il/i'Vi  = 0 is  obtained  by  adding  the  second  of  Equation  (4.3.2)  to 
the  second  of  Equation  (4.3.3),  there  being  no  contribution  from  Equation  (4.3.4).  The  full  set 
of  equations  is  shown  below: 


' 2-080 

0-146 

-0-626 

-1-121 

-0-626 

0-146 

0 “ 

f¥il 

*0-333“ 

0-146 

0-426 

-0-353 

-0-219 

0 

0 

0 

¥2 

0-194 

-0-626 

-0-353 

1-566 

-0-348 

-0-141 

0 

-0  095 

¥3 

0-333 

— 1-121 

-0-219 

-0-348 

2-863 

-0-348 

—0-219 

-0-606 

¥4 

= 

0-416 

-0-626 

0 

-0-141 

-0-348 

1-566 

-0-353 

-0-095 

¥5 

0-333 

0-146 

0 

0 

-0-219 

-0-353 

0-426 

0 

¥6 

0-194 

0 

0 

-0-095 

-0-606 

-0-095 

0 

0-191 _ 

_¥?_ 

_0- 1 94_ 

(4.4.2) 

At  this  stage  the  boundary  condition  (3.1.2)  is  applied.  Because  of  the  extreme  simplicity 
of  the  present  example  all  nodes,  save  node  4,  are  boundary  nodes  and  for  each  of  these 
¥*  = 0;  also  the  corresponding  equations  3//5¥*  = 0 are  discarded.  Hence  (4.4.2)  reduce  to 
the  single  equation 

2-863¥„  = 0-416  (4.4.3) 

with  the  solution  ¥4  = 0145.  (Incidentally,  an  analytical  solution  for  the  torsion  function  is 
given  in  Reference  6 and,  from  this,  the  exact  value  ¥4  = 0116  can  be  found.) 

Since,  over  each  element,  <£(f,  17)  is  given  by  Equation  (2.2.2)  and  since,  from  Table  4,  global 
node  4 corresponds  to  local  node  I in  all  three  elements  (coincidentally),  it  follows  that  in  each 
element,  on  setting  ¥1  = 0- 145  with  all  other  ¥<  = 0, 

l)  = 0'l 45/Mf,  ,)  = 0 • 145(1  + fXl  + i»)/4  (4.4.4) 

The  torsion  constant  Ki  for  each  element  can  be  obtained  from  Equation  (3.4.10).  For 
example 

Ka  = bi'Vi  + bn'Vn  + biuVni  + (>iv¥iV  (4.4.5) 

which,  on  extracting  the  value  for  bi  from  Equation  (4.3.1),  and  bearing  in  mind  that  only  ¥1 
is  non-zero,  reduces  to 


A„  = 01388  x 0145  = 0 0202  (4.4.6) 

The  values  of  Kb  and  Kc  turn  out  to  be  the  same,  so  that  the  torsion  constant  for  the  complete 
section  obtained  by  carrying  out  the  summation  (3.3.5)  is 

K = 0 0606  (4.4.7) 

(The  exact  answer  as  given  by  Roark8  is  0- 1044.) 


8.  Roark,  R.  J,,  Formulas  for  Stress  and  Strain,  3rd  ed.,  p.  179,  McGraw-Hill,  New  York, 
1954. 
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Finally,  the  stress,  rx,  at  node  1 will  be  calculated,  treating  node  1 as  belonging  to  element 
a.  The  general  formula  is  given  by  (3.4.11).  Since  global  node  1 corresponds  to  local  node  II 
in  element  a,  for  the  present  calculations  f = — 1,7?=  1 From  Equation  (4.2.2)  the  derivatives 
of  the  mapping  function  at  this  point  have  the  following  values: 

atj  = 0 jc,  = -0-500  (4.4.8) 

= 0 ■ 1 66  = 0 

with  the  Jacobian  J = 0 0833,  from  Equation  (2.3.3).  Also,  from  Equation  (4.4.4) 

0£(_1,  1)  = 0-145(1  + ij)/4  = 0 0727  (4.4.9) 

f,(-l,  1)  = 0 145(1  + f)/4  = 0 

Using  these  values  in  Equation  (3.4.  II),  with  K having  the  value  0 0606  from  Equation  (4.4.7), 
gives  the  result 

Tx/r=  7-19  (4.4.10) 

(The  exact  answer  as  given  in  Reference  8 is  6-38.) 

5.  ILLUSTRATIVE  EXAMPLE  USING  EIGHT-NODE  ELEMENTS— TORSION  OF 
BAR  WHOSE  SECTION  IS  A QUADRANT  OF  A CIRCLE 

5.1  General 

As  an  example  of  the  use  of  eight-node  isoparametric  elements,  the  torsion  of  a bar  whose 
cross-section  comprises  a quadrant  of  a circle  of  unit  radius  will  be  considered  (Fig,  10).  The 
finite  element  subdivision  adopted  is  shown  in  Figure  11,  the  elements  again  being  identified  as 
“a”,  “b”  and  “c”.  Also  shown  in  Figure  1 1 are  the  global  node  numbers  (1-16)  and  the  local 
node  numbers  (I— VI II)  for  each  element.  These  last  must  follow  the  sequence  indicated  in  the 
counter-clockwise  sense,  in  each  element.  The  co-ordinates  of  the  nodes  are  as  shown  in  Table 
14  below.  (As  indicated  earlier  in  Section  2.3.2  all  the  intermediate  nodes  have  been  located 
at  the  half-way  points  of  the  sides.) 


The  calculations  for  the  present  case  follow  a very  similar  pattern  to  those  described  at 
length  in  Section  4 and,  again,  all  integrations  will  be  carried  out  using  the  formula  (4.1.2). 
The  description  here  will  be  somewhat  more  concise. 
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VIII  Q 8 


Global  nodes  numbered  in  ordinary  numerals 
Local  (element)  nodes  numbered  in  Roman  numerals 


FIG.  11  FINITE  ELEMENT  SUBDIVISION,  WITH  NODAL  NUMBERS, 
FOR  ILLUSTRATIVE  EXAMPLE  OF  SECTION  5 
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5.2  Manping  Functions 

The  mapping  functions  for  the  elements  are  obtained  by  substituting  the  values  of  and 
from  Table  14  into  Equations  (2.3.12).  However,  since  these  are  not  required  in  the  calculation, 
they  will  not  be  displayed  here.  (In  any  case,  the  accuracy  of  the  mapping  for  the  curved  boundary 
has  already  been  examined  in  Section  (2.3.2).  The  derivatives  of  the  mapping  functions,  which 
are  required  at  the  integration  points  (±  l/v'3,  ± l/v'3)  are  shown  below  for  each  element; 
these  are  obtained  by  substituting  from  Table  14  into  Equations  (2.3.21).  Also  shown  are  the 
quantities  J,  A,  B , and  C for  each  element  as  calculated  from  Equations  (2.3.3)  and  (3.4.3). 


Element  a 


TABLE  15 

Derivatives  of  Mapping  Function  for  Element  a 


-l/v'3,  -l/v'3 

1/V3,  -l/v'3 

-l/v'3,  l/v'3 

l/v'3,  1/V3 

xi 

0-25 

0-25 

0-25 

0-25 

x, 

0 

0 

0 

0 

yi 

0 

0 

0 

0 

y<i 

0-25 

0-25 

0-25 

0-25 

TABLE  16 

Values  of  J,  A,  B,  C for  Element  a 


V 

1 v 3.  l/v'3 

l/v'3.  -l/v'3 

-l/v'3,  1/V3 

1/V3,  l/v'3 

J 

0 0625 

0 0625 

0 0625 

0 0625 

A 

0 0625 

0 0625 

0 0625 

0 0625 

B 

0 

0 

0 

0 

C 

0 0625 

0 0625 

0 0625 

0 0625 

The  simple  form  of  these  results  is,  of  course,  due  to  the  fact  that  element  a is  a square  so 
that  here  the  mapping  is  simply  of  one  square  onto  another. 


Element  b 


TABLE  17 

Derivatives  of  Mapping  Function  for  Element  b 


(,V  -l/v'3.  -l/v'3  l/v'3.  —l/v'3 

— l/v'3,  l/v'3 

l/v'3,  l/v'3 

x(  0-2424  0-2424 

0-1579 

0-1579 

x,  -0-0137  -0-0514 

-0-0480 

-0-1795 

yc  0-0315  0-0315 

0-0913 

0-0913 

y,  0-2789  0-3581 

0-2647 

0-3051 

TABLE  18 

Values  of  J,  A,  B,  C for  Element  b 


1/V3,  — 1/V3 


1/V3,  1/V3 


0 0462 
0 0724 
00166 
0 0332 


1/V3,  1/V3 


-0  0004 
0 0332 


Element  c 


There  is  no  need  to  do  any  calculations  for  element  c since  it  is  symmetric  with  respect  to 
element  b. 


5.3  Element  Matrices 

The  values  of  the  quantities  m,  v<,  and  w(  at  the  integration  points,  as  computed  from 
Table  3,  are  shown  in  Tables  19,  20  and  21  below. 


TABLE  19 

Values  of  U{  at  Integration  Points 


-1/V3,  -1/V3  1/V3,  -1/V3  —I/a/3,  1/a/3 


-01830 
-0  0610 
-0-6830 
-0-2276 
0-2440 
-0-3333 
0-9106 
0-3333 


0-0610 

0-1830 

0-2276 

0-6830 

-0-2440 

-0-3333 

-0-9106 

0-3333 


-0-0610 

-0-1830 

0-9106 

-0-3333 

0-2440 

0-3333 


1/V3,  l/y/3 


0-6830 

0-2276 

0-1830 

0-0610 

-0-9106 

-0-3333 

-0-2440 

0-3333 


TABLE  20 

Values  of  v<  at  Integration  Points 


-1/V3,  -1/V3  1/V3,  -I/a/3  -1/a/3,  1/V3 


1/a/3,  l/v/3 


■1830 

•2276 

•6830 

0-0610 

0-3333 

0-9106 

0-3333 

0-2440 


0-2440 

-0-3333 

0-9106 


0610 
•6830 
•2276 
•1830 
•3333 
-0-9106 
-0-3333 
-0-2440 


TABLE  21 

Values  of  w<  at  Integration  Points 


1 

-1/V3,  -1/V3 

1/V3,  -1/V3 

1 

-l/v/3,  l/v/3 

1/V3,  1/V3 

I 

Bn 

0-0962 

II 

mm 

-0-1666 

III 

WBSm 

HSSSfl 

Bfpi 

-0-0962 

IV 

BUSS 

-0-1666 

Mm 

Kssa 

01408 

0-5257 

EH 

01408 

0-5257 

0-1408 

0-5257 

0-1408 

0-1408 

H 

0 1408 

0-5257 

0-1408 

0-5257 

The  quantities  atj  and  bt  can  now  be  calculated  for  each  element  using  Equations  (3.4.7) 
and  (3.4.8)  in  conjunction  with  the  numerical  values  of  Tables  16,  18,  19,  20  and  21.  (The 
integration  formula  (4.1.2)  is,  of  course,  used.)  Now  the  ay  comprise  the  elements  of  a symmetric 
8x8  matrix.  The  results  for  each  element,  after  the  conversion  from  element  nodes  to  global 
nodes  are  shown  in  Equations  (5.3.1),  (5.3.2)  and  (5.3.3). 
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5.4  Global  Equations  and  Solution 

The  global  equations  as  given  by  Equations  (3.3.3)  are  obtained  by  the  appropriate  addition 
of  Equations  (5.3.1),  (5.3.2)  and  (5.3.3).  (For  example,  the  equation  il/i'Yn  = 0 is  obtained  by 
adding  the  first  of  (5.3.1),  the  second  of  (5.3.2),  and  the  second  of  (5.3.3).)  This  leads  to  16 
equations  in  16  unknowns.  In  the  interests  of  brevity  these  will  not  be  written  down.  After  the 
application  of  the  boundary  conditions,  which  involves  setting  *E*  = 0 for  all  k save  7,  10,  11 
and  12,  and  the  discarding  of  the  corresponding  equations  i//<>'E*  = 0,  the  reduced  set  of 
equations  for  solution  becomes 


" 4-130  0 -1-445  0-528' 

-r7  - 

' 0-330" 

0 4-130  -1-445  0-528 

0-330 

-1-445  -1-445  4-358  -3-054 

*En 

-0146 

0-528  0-528  -3  054  6-274_ 

.’Em. 

0-321. 

The  solution  of  (5.4.1)  is 

'¥^  = 0100,  'Em  = 0100,  = 0 086,  Tn  = 0 076  (5.4.2) 

Hence,  on  recalling  the  correspondence  between  global  nodes  and  element  nodes  given  in 
Table  14,  it  follows  that,  in  element  a, 

<Kf,  V)  = 0-086<?i(f  ,1 ,)  + 0-  100{?v(f,  v)  + 0 100Qvm((,  v)  (5.4.3) 

in  element  b, 

<K£,  v)  = 0 086Qn(f,  v)  + 0 076 Qv({,  r,)  + 0-  100(?vi(f,  v)  (5.4.4) 

and  in  element  c, 

<Kf,  n)  = 0 086<2iv(£,  v)  + 0-  100<?vii(f,  v)  + OO760vin(f,  i>)  (5.4.5) 

The  torsion  constant  for  each  element  is  obtained  using  Equation  (3.4.10).  For  example,  for 
element  a, 

Ka  = 0-086hi  + 0-  1006v  + 0-  lOOAvm  (5.4.6) 

and,  extracting  the  relevant  values  of  bi  from  Equation  (5.3.1),  gives 

Ka  = 0 0298 

Analogous  calculations  for  the  other  elements  give 

Kb  = Ke  = 0 0241 

Summing  the  values  for  the  elements  gives  the  torsion  constant  for  the  complete  section;  the 
result  is 


K=  0 0781  (5.4.7) 

(The  exact  answer  as  given  in  Reference  8 is  0 0825.)  The  stresses  rz  and  r„  at  node  13  will  be 
calculated,  treating  this  node  as  belonging  to  element  b.  Since  global  node  13  corresponds  to 
local  node  I in  element  b,  the  following  calculations  are  made  with  f = -q  = 1 in  Equations 
(3.4.1 1).  The  derivatives  of  the  mapping  function,  as  obtained  from  Equations  (2.3.21)  take  the 
values 

x(  = 0 1035  x,  = -0-2870  (5.4.8) 

>-{  = 01035  ^ = 0-2952 

with  the  Jacobian  J = 0 0603. 

From  Equation  (3.4.4)  and  Table  3,  it  follows  that 

= 0-086(1  + v)(2(  - q)/4  + 0-076(-fKl  + v)  - 0- 100(1  - q*)/2 

Hence, 


it( I,  D=  -0-109 


(5.4.9) 
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An  analogous  calculation  gives 


0,(1,  1)  = 0 (5.4.10) 

Substituting  from  (5.4.8),  (5.4.9)  and  (5.4.10)  into  Equations  (3.4.11),  with  K having  the  value 
0 0781  leads  to  the  result 

rx/T=-  6-67  ry/T=  +6-86  (5.4.11) 

From  considerations  of  symmetry  it  would  be  expected  that  the  two  stresses  would  be  equal. 
Actually,  if  the  calculation  be  repeated,  now  treating  node  13  as  belonging  to  element  c,  the 
results  are  reversed,  i.e. 

tx/T=  -6-86  r„/7’=+6-67  (5.4.12) 

In  general,  there  are  discontinuities  in  the  stresses  at  the  boundaries  of  adjacent  elements  and, 
clearly,  if  mean  values  be  taken  here,  the  equality  required  by  symmetry  is  restored. 

6.  OTHER  APPLICATIONS  OF  ISOPARAMETRIC  ELEMENTS 

The  two  types  of  isoparametric  elements  described  above  can  be  used  for  a wide  variety 
of  two-dimensional  problems  in  continuum  mechanics.  As  well  as  problems  in  elasticity,  prob- 
lems in  heat  conduction,  fluid  mechanics,  etc.,  can  be  solved.  The  discussion  of  Section  2 is 
applicable  in  all  cases  but,  naturally,  that  of  Section  3 must  be  replaced  by  the  appropriate 
physical  formulation. 

It  is  possible  to  extend  the  concept  to  three-dimensional  elements.  One  can  develop  an 
eight-node  plane-sided  cuboid  which  is  the  three-dimensional  analogue  of  Figure  4 and  a 
twenty-node  curved-sided  cuboid  which  is  the  analogue  of  Figure  5.  However,  particularly  in 
the  latter  case,  the  analysis  becomes  formidable. 

Reverting  again  to  the  two-dimensional  situation,  the  application  of  the  eight-node  element 
has  received  considerable  attention  for  the  determination  of  the  stress  intensity  factor  at  the  tip 
of  a crack  in  an  elastic  sheet.9'10  The  procedures  of  References  9 and  10  differ  somewhat  and  the 
following  discussion  is  based  on  the  latter  reference.  A wedge-shaped  element  (Fig.  12)  is  used 
and  this  is  obtained  by  collapsing  one  side  of  an  originally  four-sided  region  into  a single  point, 
which  is  located  at  the  crack  tip.  An  important  detail  is  that  the  intermediate  nodes  V and  VII 
must  be  placed  at  one-quarter  of  the  distance  along  each  side  from  the  crack  tip.  If  the  vertex 
angle  be  denoted  by  a,  then  the  co-ordinates  of  the  nodes  are  typically  as  shown  in  Table  22. 


TABLE  22 

Nodal  Co-ordinates  for  Crack-Tip  Element 


Node 

Xi 

y< 

I 

COS  a 

sin  a 

II 

0 

0 

III 

0 

0 

IV 

l 

0 

V 

(cos  a)/4 

(sin  a)/4 

VI 

0 

0 

VII 

1/4 

0 

VIII 

(1  + COS  a)/2 

(sin  a)/2 

9.  Henshell,  R.  D.,  and  Shaw,  K.  G.,  Crack  Tip  Finite  Elements  are  Unnecessary,  Inter.  J. 
Numerical  Methods  in  Engineering,  vol.  9,  pp.  495-507,  1975. 

10.  Hussain,  M.  A.,  Lorenson,  W.  E.,  and  Pflegl,  G.,  The  Quarter-Point  Quadratic  Element  as  a 
Singular  Element  for  Crack  Problems,  NASTRAN:  Users’  Experiences,  Fifth  Colloquium, 
NASA  TM  X-3428,  pp.  419-38,  October  1976. 
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The  region  of  Figure  12  is  mapped  on  to  that  of  Figure  3 using  the  standard  transformation 
(2.3.12).  On  substituting  the  values  of  and  yt  from  Table  22  into  (2.3.12)  it  is  found,  after 
some  trigonometric  manipulations,  that 

* = (1  + £)2{cos2(a/2)  - r,  sin2(a/2)}/4  (6.1) 

y = 0 + f)2(l  + ^){sin  (a/2)  cos  (a/2)}/4 
The  Jacobian  of  this  transformation  is  given  by 

J = (1  + f)3  sin  a/16  (6.2) 

which,  naturally,  vanishes  all  along  f — — I.  Introducing  polar  co-ordinates  r,  6 given  by 

x = r cos  6,  y = r sin  6 (6.3) 

it  is  possible  to  invert  Equations  (6.1)  to  obtain 

£ = {2 r*  cos*  ( 6 - a/2)}/{cos*  (a/2)}  - 1 (6.4) 

t)  = {tan  ( 6 — a/2)}/tan  a/2 

A displacement  component,  u,  is  taken  to  be  given  by  the  general  formula  analogous  to  Equation 
(2.2.5),  i.e. 

VIII 

« = £ £?<(£  ’))«<  (6.5) 

I = I 

From  this  last,  the  formula  for  u(f,  -q)  can  be  obtained  by  setting  un  — uni  = uvi  = 0. 
If,  in  this  formula — which  will  not  be  written  down  at  length  here — the  substitutions  (6.4)  are 
made,  it  will  be  found  that  in  the  vicinity  of  the  crack  tip,  where  r is  small,  u is  proportional 
to  r*.  This  is  the  characteristic  behaviour  required  of  a crack-tip  element.  For  further  details, 
reference  should  be  made  to  the  papers  cited  previously. 


II 
VI 

III 


I 


FIG.  12  ISOPARAMETRIC  ELEMENT  FOR  DETERMINING  STRESS  INTENSITY  FACTOR 


7.  CONCLUSIONS 

The  main  point  of  the  present  work  has  been  to  exemplify  ,by  means  of  particular  problems 
treated  in  some  detail,  the  use  of  isoparametric  finite  elements.  Such  elements  can  be  gainfully 
employed  in  a wide  variety  of  applications. 
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